% =========================================================================
% InitParameters.m - (修改版) 初始化所有系统参数和NMPC控制器
%
% 核心修改:
% 1. 扩展状态(x)和控制(u)向量的定义，以包含温度和底层流量。
% 2. 增加了详细的燃料电池物理模型参数(FCParams)，并标注了来源或假设。
% 3. 重构了NMPC控制器(nlobj)的尺寸、输出函数和约束，以适应新模型。
% =========================================================================
clear; clc; close all;

%% 1. 仿真参数
SimParams.Ts = 0.1;       % 系统采样时间 (s)
SimParams.T_sim = 50;     % 总仿真时间 (s)

%% 2. 电池模型参数 (BattParams) - 保持不变
BattParams.Capacity_As = 1 * 360; % 电池总容量 (单位: 安秒 A*s), 1Ah
BattParams.OCV_SOC_map = [
    0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7,  0.8,  0.9,  1.0;
    41.0, 41.8, 42.5, 43.2, 44.4, 45.6, 46.8, 48.0, 49.2, 50.4
];

%% 3. 燃料电池模型参数 (FCParams) - 【重大更新 + 参数缩放】
% 衰减模型系数 (保持不变)
FCParams.C1 = 1e-8;      % 负载大小相关的衰减系数 (调整后)
FCParams.C2 = 5e-7;      % 功率波动相关的衰减系数 (调整后)

% --- 新增：详细物理模型参数 (参考论文并为无人机缩放) ---

% a) 首先定义无人机级别的功率范围 (这是您之前缺少的)
FCParams.P_fc_min = 50;        % 电堆最小输出功率 (W) 【根据您的应用设定】
FCParams.P_fc_max = 500;       % 电堆最大输出功率 (W) 【根据您的应用设定】

% b) 电化学相关参数 (进行等比例缩放)
FCParams.u_H2_util = 0.95;     % 氢气利用率 (典型值) 【假设值】
FCParams.F_faraday = 96485;    % 法拉第常数 (C/mol) (物理常数)
FCParams.R_gas = 8.314;        % 理想气体常数 (J/(mol*K)) (物理常数)
FCParams.B_const = 0.05;       % 浓差极化常数 (V) (论文)
FCParams.J_max = 1.5 / (1e-2)^2; % 最大电流密度 (A/m^2), 论文中1.5 A/cm^2 (保持不变)

% --- 参数缩放的核心 ---
% 目标: 设计一个额定电压约48V的无人机电堆
target_voltage = 48; % V
rated_cell_voltage = 0.6; % V, 额定功率下的单电池电压 (经验值)
FCParams.N_cells = round(target_voltage / rated_cell_voltage); % ==> 约80个电芯

% 根据最大功率和电压，计算最大电流
max_I = FCParams.P_fc_max / target_voltage; % A
% 根据最大电流和最大电流密度，反推所需的电极面积
FCParams.A_cell_area = max_I / FCParams.J_max; % ==> 单位: m^2

% c) 压缩机相关参数 (进行缩放)
FCParams.kappa_air = 1.4;      % 空气绝热指数 (论文)
FCParams.eta_ad = 0.75;        % 压缩机绝热效率 (论文假设范围0.75-0.85) 【假设值】
FCParams.p_s = 1.013e5;        % 入口压力, 标准大气压 (Pa) (论文)
FCParams.p_a = 1.5e5;          % 阳极入口压力 (Pa) 【假设值, 需高于饱和水蒸气压】
FCParams.p_c = 1.5e5;          % 阴极入口压力 (Pa) 【假设值】
% 重新估算压缩机功耗系数。假设最大功率时，压缩机功耗占20%
max_q_H2 = (max_I * FCParams.N_cells) / (2 * FCParams.u_H2_util * FCParams.F_faraday);
max_q_air = max_q_H2 * 2.5 / 0.21; % 假设最大功率时化学计量比为2.5
P_comp_max = FCParams.P_fc_max * 0.20;
FCParams.P_comp_coeff = P_comp_max / (max_q_air^2); % ==> (W/(mol/s)^2) 【基于缩放的估算】

% d) 热力学相关参数 (进行缩放)
% 假设热容与电堆质量(正比于电芯数量和面积)成正比
scaling_factor = (FCParams.N_cells * FCParams.A_cell_area) / (332 * 0.03); % 与论文模型的尺寸比例
FCParams.C_thermal = 8000 * scaling_factor;     % ==> 系统总热容 (J/K) 【基于缩放的估算】
FCParams.a_rad_area = (0.8*0.2) * scaling_factor; % ==> 辐射散热面积 (m^2) 【基于缩放的估算】

FCParams.delta_emissivity = 0.8; % 辐射发射率 (论文)
FCParams.sigma_boltzmann = 5.67e-8; % 玻尔兹曼常数 (W/(m^2*K^4)) (论文)
FCParams.Cp_water = 75.3;      % 冷却水比热容 (J/(mol*K)) (论文)
FCParams.T0 = 298.15;          % 环境温度 (K) (论文)
FCParams.l_membrane = 1e-4;    % 质子膜厚度 (m) 【典型值, 假设】

% e) 最小流量计算 (现在可以正确计算了)
% 假设最低功率时，电堆电压较高，例如 0.8V/cell * N_cells
min_I = FCParams.P_fc_min / (0.8 * FCParams.N_cells); 
FCParams.q_H2_min = (min_I * FCParams.N_cells) / (2 * FCParams.u_H2_util * FCParams.F_faraday);
FCParams.q_air_min = FCParams.q_H2_min * 2.5 / 0.21; % 假设最小化学计量比为2.5
FCParams.q_W_min = 0; % 最小冷却水流量
%% 4. 成本函数参数 (CostParams) - 保持不变
CostParams.p_H2 = 30; 
CostParams.k_H2 = 1.67e-8; 
CostParams.C_bat_total = 500;
CostParams.N_cycle = 1000;
CostParams.C_fc_total = 2000;
CostParams.alpha_max = 20;

%% 5. 创建并配置非线性MPC控制器 (nlmpc) - 【重大更新】

% --- 定义新的模型尺寸 ---
nx = 4; % 状态量数量 [SOC, T_stack, alpha_decay, P_fc_actual]
ny = 4; % 输出量数量 [SOC, T_stack, alpha_decay, P_fc_actual]
nu = 3; % 控制量数量 [q_H2_in, q_air, q_W]

% --- 创建NMPC对象 ---
nlobj = nlmpc(nx, ny, nu);

% --- 指定预测模型 ---
nlobj.Model.StateFcn = 'prediction_model';
nlobj.Model.NumberOfParameters = 1; % 扰动 P_load

% 定义输出函数, 输出即状态
nlobj.Model.OutputFcn = @(x,u,p) x; 

% --- 设置时间步长与时域 ---
nlobj.Ts = SimParams.Ts;
nlobj.PredictionHorizon = 15;
nlobj.ControlHorizon = 1;

% --- 指定自定义成本函数 ---
nlobj.Optimization.CustomCostFcn = 'cost_function';
nlobj.Weights.OutputVariables = [0 0 0 0];
nlobj.Weights.ManipulatedVariables = [0 0 0];
nlobj.Weights.ManipulatedVariablesRate = [0.0 0.0 0.0]; % 增加对控制量变化的惩罚

% 采用离散时间
nlobj.Model.IsContinuousTime = false;
% 缩放因子的选择通常是变量的典型工作范围或最大值
nlobj.States(1).ScaleFactor = 1;      % SOC (范围 0-1)
nlobj.States(2).ScaleFactor = 350;    % T_stack (典型值 353K)
nlobj.States(3).ScaleFactor = 10;     % alpha_decay (范围 0-20)
nlobj.States(4).ScaleFactor = FCParams.P_fc_max/5; % P_fc (最大值 500W)

nlobj.ManipulatedVariables(1).ScaleFactor = max_q_H2/5; % q_H2_in
nlobj.ManipulatedVariables(2).ScaleFactor = max_q_air/5; % q_air
nlobj.ManipulatedVariables(3).ScaleFactor = 1; % q_W (最大值设为1)
% --- 设置新的约束条件 ---
% 1. 对新的控制量 u 的约束
nlobj.ManipulatedVariables(1).Min = FCParams.q_H2_min; % q_H2_in
nlobj.ManipulatedVariables(1).Max = 0.14;              % 氢气流量上限 (mol/s) 【经验值】

nlobj.ManipulatedVariables(2).Min = FCParams.q_air_min; % q_air
nlobj.ManipulatedVariables(2).Max = 0.7;               % 空气流量上限 (mol/s) (参考论文图4)

nlobj.ManipulatedVariables(3).Min = FCParams.q_W_min;  % q_W
nlobj.ManipulatedVariables(3).Max = 10;                % 冷却水流量上限 (mol/s) (参考论文图5)

% 2. 对状态量 x 的约束
nlobj.States(1).Min = 0.2;    % 状态1 (SOC)
nlobj.States(1).Max = 0.8;

nlobj.States(2).Min = 298.15; % 状态2 (T_stack) 下限，环境温度
nlobj.States(2).Max = 363.15; % 状态2 (T_stack) 上限，90摄氏度

nlobj.States(3).Min = 0;      % 状态3 (alpha_decay)

nlobj.States(4).Min = FCParams.P_fc_min; % 状态4 (P_fc)
nlobj.States(4).Max = FCParams.P_fc_max;

% 3. 对燃料电池功率变化率的约束 (使用自定义约束)
nlobj.Optimization.CustomIneqConFcn = 'power_rate_constraint';

disp('初始化完成，NMPC控制器已创建并配置完毕。');
